версия 0.4 от 05.09.2025
Хабибуллин Ринат
Работа с моделями в ходе курса может проводиться с использованием различного ПО
python для построения аналитических моделей
Для работы можно использовать jupyter notebook или google colab или любой другой сервис для работы со скважинами. Потребуется освоить базовые библиотеки python - numpy, scipy, matplotlib, mpmath, pandas. Также будут использованы некоторые другие специализированные библиотеки.
Макросы MS Excel
Часть расчетов можно выполнить с использованием макросов MS Excel, например Unifloc VBA (https://github.com/unifloc/unifloc_vba)
Гидродинамический симулятор для построения численных моделей
Базовый вариант - тНавигатор. Версия начиная от 19. Можно будет использовать другие симуляторы - РН-КИМ, Schlumberger eclipse, OPM и другие.
Специализированное ПО для анализа ГДИС
Возможно, что то появится в ходе курса.
Уравнение фильтрации многофазного потока в пористой среде - основная математическая модель пласта используемая в инженерных приложениях. Используется не само уравнение, а разнообразные частные решения, описывающие те или иные ситуации. Наиболее широко известным примером решения уравнения фильтрации, возможно является, численное решение которое строится трехмерными гидродинамическими симуляторами (eclipse, tnavigator и так далее).
Как правило, при выводе уравнения фильтрации используются следующие соотношения:
\[ u_r=-\frac{k}{\mu}\frac{dp}{dr} \tag{1}\]
уравнение неразрывности - или закон сохранения массы, показывает что в потоке не могут образовываться разрывы сплошности: \[ \frac{1}{r}\frac{\partial\left(r\rho u_r\right)}{\partial r}=-\varphi\frac{\partial p}{\partial t} \tag{2}\]
уравнение состояния - задает зависимость плотности флюида от давления: \[ c_0=\frac{1}{\rho}\frac{\partial\rho}{\partial p} \tag{3}\]
Опираясь на эти соотношения уравнение фильтрации в радиальной форме можно привести к виду для величин СИ.
\[ \frac{\partial ^2 p }{\partial r^2} + \frac{1}{r} \frac{\partial p}{\partial r} = \frac{\varphi \mu c_t}{k} \frac{\partial p}{\partial t} \tag{4}\]
Здесь используются следующие обозначения:
Иногда бывает удобно записать уравнение используя другие единицы измерения, например практические метрические - уравнение примет вид:
\[ \frac{\partial ^2 p }{\partial r^2} + \frac{1}{r} \frac{\partial p}{\partial r} = \frac{\varphi \mu c_t}{0.00036 k} \frac{\partial p}{\partial t} \tag{5}\]
Вывод уравнения фильтрации основан на следующих предположениях:
Уравнение фильтрации основано на трех основных уравнениях:
Выпишем эти уравнения, чтобы разобраться. Ориентируемся в обозначениях и подходе на [1]. Особое внимание на данном этапе стоит обратить на размерные коэффициенты в уравнениях. С одной стороны записывая уравнения в [СИ] их можно упростить и избежать части сложностей с переводными коэффициентами. С другой стороны во многих источниках переводные коэффициенты интесивно используются и это оказывается удобно при программной реализации алгоритмов. В данном описании ориентируясь на [1] приведем перевод коэффициентов в практические метрические единицы.
В дифференциальной форме, для потока в направлении \(x\)
\[ -f \dfrac{\partial \rho q_x}{\partial x} = A \dfrac{\partial(\rho \phi)}{\partial t} \tag{6}\]
где
В перечне указаны размерности для промысловых метрических и американских промысловых единиц измерений, там где они отличаются.
Для практических метрических единиц (6) можно записать в виде \[ -0.04167 \dfrac{ \partial (\rho q_x) }{\partial x} = A \dfrac{\partial(\rho \phi)}{\partial t} \tag{7}\]
Простейшим вариантом задания зависимости между измением давления и потоком флюида является линейная зависимость - закон Дарси. При этом линейный множитель в зависимости - функция как коллектора (проницаемость), так и свойств флюида (вязкость).
Закон Дарси в направлении \(x\). \[ q_x = - \dfrac{k_x A}{f\mu} \dfrac{\partial p}{\partial x} \tag{8}\]
где
В практических метрических единицах измерения Уравнение 8 будет иметь вид \[ q_x = - \dfrac{k_x A}{115.74\mu} \dfrac{\partial p}{\partial x} \tag{9}\]
Пытаясь проверить корректность вывода переводных коэффициентов необходимо учесть что для единиц измерение СИ получим \(f=1\), а необходимые переводные коэффициенты можно записать как:
В простейшем варианте, который мы будет изучать, предполагается, что сжимаемость постоянна и мала
\[ c_t = c_{formation} + c_{fluid} = const \tag{10}\]
Подставим выражение для закона Дарси Уравнение 9 в уравнение мат баланса Уравнение 7
\[ 0.04167 \dfrac{ \partial}{\partial x} { \left[ \rho \dfrac{k_x A}{115.74\mu} \dfrac{\partial p}{\partial x} \right] } = A \dfrac{\partial(\rho \phi)}{\partial t} \tag{11}\]
Упрощая Уравнение 11 получим
\[ \frac{\partial (\rho \phi)}{\partial t} = 0.00036 k_x \dfrac{ \partial}{\partial x}\left[ \dfrac{\rho}{\mu}\dfrac{\partial p}{\partial x} \right] \tag{12}\]
Уравнение 12 подходит для потока газа. Для жидкости, предполагая малую сжимаемость флюида можно его упростить до вида
\[ \frac{\partial p}{\partial t} = 0.00036 \dfrac{k_x}{\phi \mu c_t} \dfrac{ \partial^2 p}{\partial x^2} \tag{13}\]
Напомним, здесь
Приведенное Уравнение 13 соответствует линейному потоку в пласте. Практический интерес представляет решение для радиального потока
Закон сохранения массы или принцип неразрывности можно выразить в радиальной форме следующим соотношением
\[ \frac{1}{r}\frac{\partial\left(r\rho u_r\right)}{\partial r}=-\varphi\frac{\partial p}{\partial t} \tag{14}\]
Принцип неразрывности показывает, что для определенного объема пласта, масса флюида которое втекла в контрольный объем пласта минус масса которая вытекла равна массе которая накопилась в объеме.
Сохранение импульса или закон Дарси можно выразить соотношением
\[ u_r=-\frac{k}{\mu}\frac{d^p}{dr} \tag{15}\]
Закон Дарси здесь используется как псевдоустановившаяся аппроксимация обобщенного уравнения сохранения импульса (то есть слагаемым отвечающим за накопление импульса пренебрегаем). Это предположение справедливо если не учитывать возмущения давления в среде двигающиеся со скоростью звука. Все изменения давления описываемые моделью связаны с локальными изменениями градиента давления за счет ламинарного режима потока. Хотя далее в модели будет учтена сжимаемость системы всеми “звуковыми” эффектами в системе мы пренебрегаем.
Поток предполагается горизонтальным, поэтому давление \(p\) может быть использовано в качестве потенциала потока, гравитационными силами пренебрегаем.
Комбинируя Уравнение 14 , Уравнение 15 получим:
\[ \frac{1}{r}\frac{\partial\left( \dfrac{r\rho k}{\mu}\dfrac{\partial p}{dr}\right)}{\partial r}=\varphi\frac{\partial \rho}{\partial t} \tag{16}\]
Уравнение 16 – дифференциальное уравнение в частных производных описывающее нестационарный поток однофазного флюида в пористой среде при ламинарном потоке.
Вообще говоря приведенное уравнение является нелинейным, так как плотность \(\rho = \rho(p)\) и вязкость \(\mu = \mu(p)\) являются функциями давления. Уравнение содержит две зависимые переменные - давление \(p\) и плотность \(\rho\). Поэтому для его решения необходимо задать еще одно соотношение, каковым может быть уравнение состояния флюида связывающее плотность флюида и давление \(\rho = \rho(p)\).
Нестационарное поведение давления в пласте связано со сжимаемостью системы. При изменении давления в какой то точке, часть флюида сжимается, происходит накопление или отдача флюида, что вызывает задержку в распространении изменения флюида. Несмотря на то, что сжимаемость флюидов и породы малы и во многих случаях ими можно пренебречь, это не верно для пластовых систем для добычи нефти. Большие объемы пласта и флюидов и высокие давления компенсируют малость сжимаемости и требуют ее учета.
Для однофазного флюида разумным является предположение постоянства сжимаемости.
Сжимаемость можно определить как
\[ c=-\frac{1}{V} \left( \frac{ \partial V}{ \partial p} \right) \]
учтем, что
\[ \rho = \frac{m}{V} \]
тогда получим
\[ c=\frac{1}{\rho} \left( \frac{ \partial \rho}{ \partial p} \right) \]
Для флюида с постоянной сжимаемостью, проинтегрировав приведенное уравнение можно получить
\[ \rho = \rho_i e^{c(p-p_i)} \]
где \(\rho_i\) плотность флюида при некотором заданном давлении \(p_i\)
Продифференцировав выражение для плотности по времени получим
\[ c \rho \frac{\partial p}{\partial t} = \frac{\partial \rho}{\partial t} \]
Подставив это выражение в ранее полученное уравнение фильтрации (Уравнение 16) получим
\[ \frac{1}{r}\frac{\partial\left( \dfrac{r\rho k}{\mu}\dfrac{\partial p}{dr}\right)}{\partial r}=\varphi c \rho \frac{\partial p}{\partial t} \tag{17}\]
Приведенное дифференциальное уравнение в частных производных все еще нелинейно, поскольку зависит от плотности \(\rho\)
Если пористость не является постоянной величиной и меняется с давлением, тогда слагаемое отвечающее за накопление флюида в пласте можно выразить как
\[ \frac{\partial \varphi \rho}{\partial t} = \varphi \frac{\partial \rho}{\partial t}+ \rho \frac{\partial \varphi }{\partial t} = \varphi c_l \rho \frac{\partial p}{\partial t} + \rho \frac{\partial \varphi }{\partial t} \]
где \(c_l\) сжимаемость жидкости.
Определим сжимаемость породы как
\[ c_f = \frac{1}{\varphi} \frac{\partial \varphi}{\partial p} \]
тогда
\[ \frac{\partial \varphi \rho}{\partial t} = \varphi \rho (c_l + c_f) \frac{\partial p}{\partial t} \]
хотя пористость здесь является функцией давления - в первом приближении мы можем считать ее константой равной пористости при некотором среднем давлении в пласте. Это справедливо для маленькой сжимаемости породы, что верно почти всегда.
Уравнение для сжимаемости можно еще уточнить, учтя что в пласте могут находится различные флюиды - вода и нефть с насыщенностями \(s_w\) и \(s_o\)
тогда
\[ c_l = s_o c_o + s_w c_w \]
тогда можно ввести общую сжимаемость системы
\[ c_t = c_l + c_f = s_o c_o + s_w c_w + c_f \]
Заметим, что проницаемость k в законе Дарси это не абсолютная проницаемость, но относительная проницаемость по нефти при насыщенности водой соответствующей связанной воде.
\[ k= k_o (s_{wc}) \]
Раскрыв производную в левой части Уравнение 17 и предположив, что \(\dfrac{\partial p}{\partial r}\) мало а следовательно слагаемым \(r \rho c_t \left( \dfrac{\partial p}{\partial r} \right)^2\) можно пренебречь, получим линеаризованное уравнения фильтрации
\[ \frac{\partial^2 p}{\partial r^2} + \frac{1}{r} \frac{\partial p}{\partial r}= \frac{\varphi \mu c_t}{k} \frac{\partial p}{\partial t} \tag{18}\]
Часто для анализа уравнений неустановившейся фильтрации используются безразмерные переменные. Мы будем использовать переменные в виде:
\[ r_D = \frac{r}{r_w} \]
\[ t_D = \frac{kt}{\varphi \mu c_t r_w^2} \]
\[ p_D = \frac{2 \pi kh}{q_s B \mu} \left( p_i - p_{wf} \right) \]
\[ q_D = \frac{q}{q_{ref}} \]
Здесь использованы единицы измерения СИ.
Использование безразмерных переменных позволяет упростить уравнение фильтрации, которое примет вид
\[ \frac{\partial p_D}{ \partial t_D} = \frac{1}{r_D} \frac{ \partial{ \left( r_D \dfrac{\partial p_D}{ \partial r_D} \right) } }{ \partial{r_D} } \]
\[ \frac{\partial p_D}{ \partial t_D} = \dfrac{1}{r_D} \left[ \dfrac{\partial}{\partial r_D} \left( r_D \dfrac{ \partial p_D} {\partial r_D} \right) \right] \]
Решение этого уравнения - функция безразмерного давления от безразмерных времени и расстояния \(p_D(r_D, t_D)\)
Для практических расчетов удобнее бывает использовать безразмерные переменные полученные для практических метрических единиц измерения. \[ r_D = \frac{r}{r_w} \]
\[ t_D = \frac{0.00036 kt}{\varphi \mu c_t r_w^2} \]
\[ p_D = \frac{kh}{ 18.4 q_s B \mu} \left( p_i - p_{wf} \right) \]
\[ q_D = \frac{q}{q_{ref}} \]
Здесь использованы практические метрические единицы измерения.
Уравнение фильтрации для радиального потока в линеаризованном виде можно записать в виде
\[ \frac{\partial p}{\partial t} = 0.00036 \dfrac{k}{\phi \mu c_t} \dfrac{1}{r} \left[ \dfrac{\partial}{\partial r} \left( r \dfrac{ \partial p} {\partial r} \right) \right] \]
Использование безразмерных переменных позволяет упростить уравнение фильтрации, которое примет вид
\[ \frac{\partial p_D}{ \partial t_D} = \dfrac{1}{r_D} \left[ \dfrac{\partial}{\partial r_D} \left( r_D \dfrac{ \partial p_D} {\partial r_D} \right) \right] \]
Решение этого уравнения - функция безразмерного давления от безразмерных времени и расстояния \(p_D(r_D, t_D)\)
Несмотря на простоту определений безразмерных переменных их часто приходится применять при проведении расчетов. Поэтому в надстройке Unifloc VBA реализован набор функций расчета безразмерных переменных.
Эти функции начинаются с префикса transient_def.
transient_def_cd
transient_def_cs_1atm
transient_def_td
transient_def_t_day
transient_def_pd
transient_def_pwf_atmaОписания функций и из аргументов можно найти в руководстве пользователя Unifloc VBA
Широкое распространение на практике получили стационарные решения уравнения фильтрации. Приведем некоторые из них.
Рассматривается самая простая модель работы добывающей скважины - радиальная стационарная фильтрация в однородном изотропном пласте круговой формы. Скважина находится в центре пласта (Рисунок 1). На границе пласта поддерживается постоянное давление. Фактически это означает, что через границу пласта идет поток жидкости, уравновешивающий дебит скважины.
Решение можно получить как решения уравнения фильтрации, учитывая стационарность потока
\[ \frac{\partial ^2 p }{\partial r^2} + \frac{1}{r} \frac{\partial p}{\partial r} = 0 \tag{19}\]
Или можно получить его непосредственно из закона Дарси, который должен быть приведен к радиальной форме и в таком варианте известен как Формула Дюпюи. Для приведенной конфигурации можно записать закон Дарси в форме.
\[ u_r=\frac{q}{2\pi rh}=\frac{k}{\mu}\frac{dP}{dr} \]
Проинтегрировав выражение по замкнутому контуру радиуса \(r_e\) вокруг скважины получим выражение известное как формула Дюпюи
\[ q=\frac{2\pi kh\left(P_e-P_w\right)}{\mu\left(\ln{\dfrac{r_e}{r_w}}\right)} \tag{20}\]
В приведенном выражения использованы единицы СИ.
Здесь
На практике часто бывает удобнее пользоваться значениями в практических метрических единицах измерения.
\[ q=\frac{kh\left(P_e-P_w\right)}{ 18.4 \mu\left(\ln{\dfrac{r_e}{r_w}}\right)} \tag{21}\]
где
Далее если не указано особо будем использовать практические метрические единицы.
Скин-фактор — гидродинамический параметр, характеризующий дополнительное фильтрационное сопротивление течению флюидов в околоскважинной зоне пласта, приводящее к изменению добычи (дебита) по сравнению с совершенной (идеальной) скважиной. Скин-фактор может приводить как к снижению дебита (например при загрязнении ПЗС), так и увеличению (образование высокопроводящих каналов в ПЗС).
Концепция скин-фактора получила широкое распространение на практике. Все инженеры-нефтяники знают этот параметр и оперируют им на практике.
Изначально скин-фактор был введен как параметр учитывающий изменение проницаемости (загрязнение) призабойной зоны при расчете производительности скважины. Такое загрязнение может быть вызвано различными причинами:
Для модели загрязненной призабойной зоны величину скин-фактора можно выразить формулой Хокинса . Скин-фактор для плоскорадиального установившегося потока несжимаемой жидкости:
\[ S =\left( \frac{k}{k_s} -1\right)\ ln\frac{r_s}{r_w} \tag{22}\]
здесь:
Концепция скин-фактора оказалась удобной для описания характеристики соединения скважины и пласта и была распространена на другие случаи, когда производительность скважины могла отличаться от производительности идеальной скважины:
Для многих подобных случаев предположение о радиальном притоке к скважине не верно, но величину скин-фактора используют, так как она позволяет сравнить производительность скважины со сложным заканчиванием с простой вертикальной скважиной. В таких случая говорят о псевдорадиальном скин-факторе - такой величине скин-фактора \(S\), которая обеспечила бы такую же производительность для вертикальной скважины полностью вскрывающей пласт.
Для стационарной радиальной модели притока учет скин-фактора приведен к следующим соотношениям:
\[ (P_e - P_{wf}) = \frac{18.41\mu q }{\ k h}(\ln\frac{r_e}{r_w}+S) \tag{23}\]
\[ q=\frac{kh\left(P_e-P_w\right)}{ 18.41 \mu\left(\ln{\dfrac{r_e}{r_w}} + S\right)} \tag{24}\]
В приведенное выражение входит значение давления на контуре, которым не всегда бывает удобно пользоваться. В практических случаях значение на контуре трудно оценить, контур зоны дренирования может значительно отличаться от кругового, да и радиус оценить может быть сложно. Удобнее пользоваться средним давлением в зоне дренирования \(\bar{P}\), которое может быть оценено по материальному балансу (смотри рисунок \(\ref{ris:radial_inflow_steady_state_average_pressure}\)).
Вывод уравнения фильтрации для постоянного давления на границе с использованием среднего давления
В этом случае выражение для дебита примет вид \[ q=\frac{kh\left( \bar{P}-P_w\right)}{ 18.41 \mu\left(\ln{\dfrac{r_e}{r_w}} - \dfrac{1}{2}+ S \right)} \]
Схема модели радиального притока для условия непротекания на круговой границе приведена на Рисунок 3
При условии непротекания давления на границе условия стационарности (неизменности давления) не достигаются. При работе скважины с постоянным дебитом забойное давление будет постоянно снижаться. Однако начиная с некоторого момента, когда влияние скважины достигнет границ - давление в всей области дренирования начнет снижаться равномерно (смотри Рисунок 4).
Такой режим, при котором забойное давление меняется, но перепад давления \(P_e - P_w\) остается постоянным называют псевдо-установившимся режимом работы (pss - pseudo steady state).
Для псевдо-установившегося режима можно записать выражение
\[ q=\frac{kh\left(P_e-P_w\right)}{ 18.41 \mu\left(\ln{\dfrac{r_e}{r_w}} - \dfrac{1}{2} + S\right)} \] где
Аналогично случаю для постоянного давления на границе можно переписать выражение с использованием среднего давления в области дренирования.
\[ q=\frac{kh\left( \bar{P}-P_w\right)}{ 18.41 \mu\left(\ln{\dfrac{r_e}{r_w}} - \dfrac{3}{4}+ S \right)} \]
Вывод уравнений для псевдо-установившегося режима работы
Здесь уравнения и методы расчета для горизонтальных, наклонно направленных скважин, скважин с ГРП, горизонтальных скважин с МГРП.
\[ q=\frac{kh\left( \bar{P}-P_w\right)}{ 18.41 \mu\left(\ln{\dfrac{2.2458 A}{C_A r_w^2}} + S \right)} \]
%todo -надо переделать рисунок и добавить ссылки на источники
Безразмерное время достижения псевдо-установившегося режима притока, определяемое видом резервуара
\[ t_{DApss} = \frac{kt_{pss}}{\varphi \mu c_t A} \]
Простое решение для задачи стационарного притока к вертикальной скважине в однородном изотропном пласте круговой формы с постоянным давлением на границе имеет вид
\[ Q=\dfrac{kh}{18.41\mu B} \dfrac{P_{res}-P_{wf}}{\ln \dfrac{r_e}{r_w} + S} \]
где:
Это решение известно как закон Дарси или формула Дюпюи.
Выражение можно переписать в виде
\[ P_{r} = P_{res} - 18.41\dfrac{ Q\mu B }{kh} \left[ \ln\dfrac{r_e}{r} +S \right] \]
который удобен для расчёта распределения давления в пласте \(P_r\) на произвольном расстоянии от скважины \(r\). В выражении (2) задано граничное значение давления \(p_e\) на контуре \(r_e\). Расчёт позволит найти любое значение внутри контура, в том числе и забойное давление \(P_{wf}\) на \(r=r_w\)
Выражение можно переписать \[ P_{r} = P_{wf} + 18.41\dfrac{ Q\mu B }{kh} \left[ \ln\dfrac{r}{r_w} +S \right] \]
где по известному дебиту и забойному давлению можно найти давление в пласте. При известном пластовом давлении можно оценить радиус контура на котором оно достигается.
Расчет на python может быть реализован на основе библиотек numpy scipy
"""
Импортируем библиотеки для расчетов.
numpy используем для работы с массивами
и подготовки данных для построения графиков.
matplotlib используем для построения графиков
scipy для решения линейных уравнений
"""
import numpy as np
import matplotlib.pyplot as plt
import scipyДля удобства дальнейшего изложения и использования расчетных функций при создании функций и переменных на языке python названия формируются по следующим принципам:
"""
Определим функции для расчета стационарного решения
"""
def dp_ss_atm(q_liq_sm3day = 50,
mu_cP = 1,
b_m3m3 = 1.2,
kh_mDm = 40,
r_e_m = 240,
r_m = 0.1):
"""
функция расчета перепада давления в произвольной точке пласта
на расстоянии r_m от центра скважины для стационарного решения
- q_liq_sm3day - дебит жидкости на поверхности в стандартных условиях
- mu_cP - вязкость нефти (в пластовых условиях)
- B_m3m3 - объемный коэффициент нефти
- kh_mDm - kh пласта
- r_e_m - радиус контрура питания, м
- r_m - расстояние на котором проводится расчет, м
"""
return 18.42 * q_liq_sm3day * mu_cP * b_m3m3/ kh_mDm * np.log(r_e_m/r_m)
def p_ss_atma(p_res_atma = 250,
q_liq_sm3day = 50,
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10,
r_e_m = 240,
r_m = 0.1):
"""
функция расчета давления в произвольной точке пласта
на расстоянии r_m от центра скважины для стационарного решения
- p_res_atma - пластовое давление, давление на контуре питания
- q_liq_sm3day - дебит жидкости на поверхности в стандартных условиях
- mu_cP - вязкость нефти (в пластовых условиях)
- B_m3m3 - объемный коэффициент нефти
- k_mD - проницаемость пласта
- h_m - мощность пласта, м
- r_e_m - радиус контрура питания, м
- r_m - расстояние на котором проводится расчет, м
"""
return p_res_atma - dp_ss_atm(q_liq_sm3day = q_liq_sm3day,
mu_cP = mu_cP,
b_m3m3 = b_m3m3,
kh_mDm = k_mD * h_m,
r_e_m = r_e_m,
r_m = r_m)Функции расчетов могут быть использованы для построения графиков, например с использованием matplotlib
"""
Построим график распределения давления в пласте
"""
# формируем массив расстояний для которых будем проводить расчет
r_arr = np.logspace(np.log10(0.1), np.log10(240), 500)
# рассчитываем массив давлений на соответствующих расстояниях
# для расчета используется векторный расчет numpy
# для примера показана передача всех аргументов созданной функции
p_arr = p_ss_atma(p_res_atma = 250,
q_liq_sm3day = 50,
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10,
r_e_m = 240,
r_m = r_arr)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))
# рисуем график в обычных координатах
ax1.plot(r_arr, p_arr) # команда отрисовки графика по заданным массивам
ax1.plot(-r_arr, p_arr)
# настраиваем график
ax1.set_title('Воронка депрессии - распределение давления в пласте')
ax1.set_xlabel('r, m')
ax1.set_ylabel('p, atma')
# рисуем график в логарифмических координатах
ax2.plot(r_arr, p_arr) # команда отрисовки графика по заданным массивам
ax2.plot(-r_arr, p_arr)
# настраиваем график
ax2.set_title('Воронка депрессии - распределение давления в пласте')
ax2.set_xlabel('r, m')
ax2.set_ylabel('p, atma')
ax2.set_xscale('symlog', linthresh=0.1, linscale=0.5)
plt.show()
"""
Построим карту или сетку распределения давления в пласте
с использованием векторных функций numpy
"""
# зададим параметры воронки депрессии
pres = 250
r_e = 300
# зададим координатную сетку основываясь на параметрах
x = np.linspace(-r_e, r_e, 300)
y = np.linspace(-r_e, r_e, 300)
# рассчитаем вспомогательные вектора для построения сетки
xv, yv = np.meshgrid(x, y)
# зададим координаты скважины
xwell1 = 0
ywell1 = 0
# рассчитаем значение давлений во всех точках сетки
# расчет ведется для матрицы координата с использованием векторных возможностей numpy
p_mesh = p_ss_atma(r_m=((xv-xwell1)**2 + (yv-ywell1)**2)**0.5, p_res_atma=pres, r_e_m=r_e)
# удалим значения за контуром, так как в данном случае они не имеют смысла
p_mesh[np.where(p_mesh > pres)] = pres
# построим отображение в виде контурной карты
plt.rcParams["figure.figsize"] = (7,7)
plt.contour(x, y, p_mesh, levels=100)
plt.show()
Для построения карты распределения давлений в пласте полезно вспомнить, что расстояние от скважины с координатами \((x_{well}, y_{well})\) до произвольной точки пласта с координатами \((x,y)\) можно найти по формуле \[ r=\sqrt{ (x-x_{well})^2 + (y-y_{well})^2 } \]
Тогда выражение для расчета давления в любой точке пласта примет вид
\[ P_{r} = P_{res} - 18.41\dfrac{ Q\mu B }{kh} \left[ \ln\dfrac{r_e}{\sqrt{ (x-x_{well})^2 + (y-y_{well})^2 }} +S \right] \]
Простой вариант расчета - можно создать пустую матрицу со значениями давления по сетке и перебирая все точки на сетке/матрице рассчитать давления
"""
Построим карту или сетку распределения давления в пласте
для ряда/галереи скважин
"""
# зададим параметры воронки депрессии
pres = 250
r_e = 300
# зададим координаты скважины (всего будет 5 скважин у которых меняется только х координата)
xwell = [-100, -50, 0 , 50, 100]
ywell = 0
# зададим координатную сетку основываясь на параметрах
x = np.linspace(-r_e, r_e, 300)
y = np.linspace(-r_e, r_e, 300)
# рассчитаем вспомогательные вектора для построения сетки
xv, yv = np.meshgrid(x, y)
# зададим пустой список матриц с перепадами давлений от каждой отдельной скважины
p_mesh_i=[]
# для каждой скважины найдем ее влияние на давления
for xi in xwell:
# рассчитаем матрицу расстояний от элементов сетки до скважины i
r_well = ((xv-xi)**2 + (yv-ywell)**2)**0.5
# рассчитаем значение давлений во всех точках сетки для скважины i
p_mesh_i_ = p_ss_atma(r_m=r_well, p_res_atma=pres, r_e_m=r_e)
# удалим значения за контуром, так как в данном случае они не имеют смысла
p_mesh_i_[np.where(p_mesh_i_ > pres)] = pres
# сохраним влияние скважины i в списке матриц влияния отдельных скважин
p_mesh_i.append(p_mesh_i_)
# найдем суммарное влияние все скважин
p_mesh_sum = 0
for p_mesh_i_ in p_mesh_i:
# найдем сумму решений, помним что суммировать можно депрессии
p_mesh_sum = pres - p_mesh_i_ + p_mesh_sum
p_mesh_sum = pres - p_mesh_sum
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))
# построим отображение в виде контурной карты
ax1.contourf(x, y, p_mesh_sum, levels=100)
# построим отображение в виде контурной карты
ax2.contour(x, y, p_mesh_sum, levels=15)
plt.show()
"""
Построим график изменения давления по линии расположения скважин
"""
xline = np.linspace(-300, 300, 10000)
yline = 0
pline = np.zeros_like(xline)
for xi in xwell:
r_well = ((xline-xi)**2 + (yline-ywell)**2)**0.5
r_well = np.where(r_well < 0.1, 0.1, r_well)
pline = pline + p_ss_atma(r_m=r_well,
p_res_atma=0,
r_e_m=r_e,
q_liq_sm3day = 50,
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10)
pline = pres + pline
plt.plot(xline, pline)
plt.show()
Для стационарного решения работает принцип суперпозиции - сумма двух решений также будет решением, это позволяет построить карту для нескольких скважин. Давление в любой точке пласта можно найти по формуле
\[ P_{res} - P_{x,y} = \sum_{i} 18.41\dfrac{ Q_i\mu B }{kh} \left[ \ln\frac{r_e}{\sqrt{ (x-x_{w.i})^2 + (y-y_{w.i})^2 }} +S \right] \]
Выражение справедливо только если \(\sqrt{ (x-x_{w.i})^2 + (y-y_{w.i})^2 }< r_e\).
При наличии нескольких скважин можно записать выражение для оценки забойных давлений скважин
\[ P_{res} - P_{wf.j} = \sum_{i} 18.41\dfrac{ Q_i\mu B }{kh} \left[ \ln\dfrac{r_e}{\sqrt{ (x_{w.j}-x_{w.i})^2 + (y_{w.j}-y_{w.i})^2 }} +S \right] \]
Если считать забойные давления \(P_{wf.j}\) известными а дебиты скважин \(Q_i\) не известными, тогда выражение (6) можно рассматривать как систему линейных алгебраических уравнений вида
\[ AX = B \]
Где \[ A_{[i,j]} = 18.41\dfrac{ \mu B }{kh} \left[ \ln\dfrac{r_e}{\sqrt{ (x_{w.j}-x_{w.i})^2 + (y_{w.j}-y_{w.i})^2 }} +S \right] \]
\[ B_{[j]}=P_{res} - P_{wf.j} \]
такую систему можно решить например с использованием пакета scipy.linalg
# зададим параметры воронки депрессии
pres = 250
r_e = 300
# зададим координатную сетку основываясь на параметрах
x = np.linspace(-r_e, r_e, 500)
y = np.linspace(-r_e, r_e, 500)
# рассчитаем вспомогательные вектора для построения сетки
xv, yv = np.meshgrid(x, y)
# зададим координаты скважины (всего будет 5 скважин у которых меняется только х координата)
xwell = [-100, -50, 0 , 50, 100]
ywell = [0,0,0,0,0]
pwell = [100, 100, 100, 100, 100] # забойные давления скважин
qwell = [0,0,0,0,0] # дебиты скважин, должны быть найдены
# создадим заготовки матрицы А и вектора В по формулам (4) и (5)
A = np.zeros((5,5))
B = np.zeros(5)
# сформируем матрицу А
for i,v in enumerate(xwell):
for j,v in enumerate(ywell):
# найдем расстояния от одной скважины до другой
# ищем расстояния между центрами скважин
r_ij = ((xwell[i]-xwell[j])**2 + (ywell[i]-ywell[j])**2)**0.5
# если расстояние 0 значит ищем влияние скважины саму на себя
# тогда подставляем радиус скважины
if r_ij == 0:
r_ij = 0.1
# чтобы воспользоваться ранее заданной формулой в виде функции
# вызовем ее с нулем пластовым давлением и единичным дебитом
A[i,j] = - p_ss_atma(p_res_atma = 0,
q_liq_sm3day = 1,
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10,
r_e_m = 240,
r_m = r_ij)
B[i] = pres - pwell[i]
# найдем решение = значения дебитов при которых забойные равны заданным значениям
qwell = scipy.linalg.solve(A,B)
# напечатаем найденные дебиты
print(qwell)
# зададим пустой список матриц с перепадами давлений от каждой отдельной скважины
p_mesh_i=[]
# для каждой скважины найдем ее влияние на давления
for i,xi in enumerate(xwell):
# рассчитаем матрицу расстояний от элементов сетки до скважины i
r_well = ((xv-xwell[i])**2 + (yv-ywell[i])**2)**0.5
# для красоты отрисовки карты давлений пренебрежем значениями в радиусе одного метра от скважин
r_well[r_well<1]=0.1
# рассчитаем значение давлений во всех точках сетки для скважины i
p_mesh_i_ = p_ss_atma(r_m=r_well, p_res_atma=pres, r_e_m=r_e, q_liq_sm3day=qwell[i])
# удалим значения за контуром, так как в данном случае они не имеют смысла
p_mesh_i_[np.where(p_mesh_i_ > pres)] = pres
# сохраним влияние скважины i в списке матриц влияния отдельных скважин
p_mesh_i.append(p_mesh_i_)
# найдем суммарное влияние все скважин
p_mesh_sum_new = 0
for p_mesh_i_ in p_mesh_i:
# найдем сумму решений, помним что суммировать можно депрессии
p_mesh_sum_new = pres - p_mesh_i_ + p_mesh_sum_new
p_mesh_sum_new = pres - p_mesh_sum_new
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))
# построим отображение в виде контурной карты
ax1.contourf(x, y, p_mesh_sum_new, levels=100)
# построим отображение в виде контурной карты
ax2.contour(x, y, p_mesh_sum_new, levels=15)
plt.show()[263.52443121 214.68400868 202.93807722 214.68400868 263.52443121]
"""
Построим график изменения давления по линии расположения скважин
"""
xline = np.linspace(-300, 300, 10000)
yline = 0
pline = np.zeros_like(xline)
for i,xi in enumerate(xwell):
r_well = ((xline-xi)**2 + (yline-ywell[i])**2)**0.5
r_well = np.where(r_well < 0.1, 0.1, r_well)
pline = pline + p_ss_atma(r_m=r_well,
p_res_atma=0,
r_e_m=r_e,
q_liq_sm3day=qwell[i],
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10)
pline = pres + pline
plt.plot(xline, pline)
plt.show()
Классы в python помогают задавать сложные структуры данных (инкапсуляция) и работать с ними как с едиными объектами. Скважины с их параметрами, а также группы скважин могут быть заданы с использованием классов.
Класс для расчета влияния одной скважины - объединяет все необходимые данные для скважины - позволяет рассчитать депрессию и давления в любой точке пласта с учетом радиуса влияния и радиуса скважины
class Well:
"""
класс Well описывает влияние одной скважины на поле давления
используется стационарное решение
"""
def __init__(self, name, xw_m, yw_m) -> None:
self.x_m = xw_m
self.y_m = yw_m
self.name = name
self.kh_mDm = 10.0
self.b_m3m3 = 1.0
self.mu_cP = 1.0
self.skin = 0
self.r_e_m = 250.0
self.r_w_m = 0.1
self.q_liq_sm3day = 10
def calc_dp_ss_atm(self, x_m, y_m):
"""
расчет перепада давления в произвольной точке пласта
с заданными координатами
для стационарного решения уравнения фильтрации
x_m, y_m координаты в которых рассчитывается давление
"""
r_m = np.sqrt( (self.x_m - x_m)**2 + (self.y_m - y_m)**2 )
r_m = np.where(r_m >= self.r_e_m, self.r_e_m ,r_m)
r_m = np.where(r_m < self.r_w_m, self.r_w_m ,r_m)
return dp_ss_atm(q_liq_sm3day=self.q_liq_sm3day,
mu_cP=self.mu_cP,
b_m3m3=self.b_m3m3,
kh_mDm=self.kh_mDm,
r_e_m=self.r_e_m,
r_m=r_m)
def calc_dp_ss_well_atm(self, w):
return self.calc_dp_ss_atm(w.x_m, w.y_m)
def calc_p_ss_atma(self, x_m, y_m, p_res_atma):
"""
функция расчета давления в произвольной точке пласта
с заданными координатами
для стационарного решения уравнения фильтрации
x_m, y_m координаты в которых рассчитывается давление
"""
return p_res_atma - self.calc_dp_ss_atm(x_m, y_m)
class Wells:
"""
класс для группы взамодействующих скважин
"""
def __init__(self, name:list, xw_m:list, yw_m:list,
kh_mDm=10, b_m3m3=1,
mu_cP=1, r_e_m =250, r_w_m=0.1) -> None:
self.well_dict: dict[str, Well] = {}
for (n, x, y) in zip(name, xw_m, yw_m):
self.well_dict.update({n: Well(n,x,y)})
self.set_kh_mDm(kh_mDm)
self.set_pvt(b_m3m3=b_m3m3, mu_cP=mu_cP)
self.set_r_e_m(r_e_m)
self.set_r_w_m(r_w_m)
self.p_res_atma = 250
def set_qliq_sm3day(self, q_liq_list_sm3day:list):
for (wname, q) in zip(self.well_dict, q_liq_list_sm3day):
self.well_dict[wname].q_liq_sm3day = q
def set_kh_mDm(self, kh_mDm):
for wn in self.well_dict:
self.well_dict[wn].kh_mDm = kh_mDm
def set_r_e_m(self, r_e_m):
for wn in self.well_dict:
self.well_dict[wn].r_e_m = r_e_m
def set_r_w_m(self, r_w_m):
for wn in self.well_dict:
self.well_dict[wn].r_w_m = r_w_m
def set_pvt(self, mu_cP, b_m3m3):
for wn in self.well_dict:
self.well_dict[wn].b_m3m3 = b_m3m3
self.well_dict[wn].u_cP = mu_cP
def calc_p_ss_atma(self, x_m, y_m):
dp = 0
for wn in self.well_dict:
dp = dp + self.well_dict[wn].calc_dp_ss_atm(x_m, y_m)
return self.p_res_atma - dp
def estimate_qliq_from_pwf(self, wn_list:list, pwf_list:list):
dim = len(wn_list) # размерность матрицы для расчета дебитов
# создадим заготовки матрицы А и вектора В по формулам (4) и (5)
A = np.zeros((dim,dim))
B = np.zeros(dim)
for i, wni in enumerate(wn_list):
self.well_dict[wni].q_liq_sm3day = 1
for j, wnj in enumerate(wn_list):
A[i,j] = self.well_dict[wni].calc_dp_ss_well_atm(self.well_dict[wnj])
B[i] = self.p_res_atma - pwf_list[i]
q_list = scipy.linalg.solve(A,B)
for i, wni in enumerate(wn_list):
self.well_dict[wni].q_liq_sm3day = q_list[i]
xlist = [300, -300, 50]
ylist = [0, 0, 200]
names = ['1', '2', '3']
qlist = [1, -20, 10]
ww = Wells(names, xlist, ylist, r_e_m=500)
ww.set_qliq_sm3day(qlist)
xp = np.linspace(-500, 500, 500)
yp = np.linspace(-500, 500, 500)
xv, yv = np.meshgrid(xp, yp)
p_mesh = ww.calc_p_ss_atma(xv, yv)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))
# построим отображение в виде контурной карты
ax1.contour(xp, yp, p_mesh, levels=50)
for x,y,n in zip(xlist, ylist, names):
ax1.text(x+10, y+10, n)
ax1.scatter(xlist, ylist, s =100)
ax2.streamplot(xp, yp, -np.gradient(p_mesh,axis=1),
-np.gradient(p_mesh,axis=0),
density=2)
for x,y,n in zip(xlist, ylist, names):
ax2.text(x+10, y+10, n)
ax2.scatter(xlist, ylist, s =100)
plt.show()
xlist = [300, -300, 50]
ylist = [0, 0, 200]
names = ['1', '2', '3']
qlist = [10, -20, 10]
plist = [50, 150, 50]
ww = Wells(names, xlist, ylist, r_e_m=500)
ww.estimate_qliq_from_pwf(names, plist)
xp = np.linspace(-500, 500, 500)
yp = np.linspace(-500, 500, 500)
xv, yv = np.meshgrid(xp, yp)
p_mesh = ww.calc_p_ss_atma(xv, yv)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))
# построим отображение в виде контурной карты
ax1.contour(xp, yp, p_mesh, levels=50)
for x,y,n in zip(xlist, ylist, names):
ax1.text(x+10, y+10, n)
ax1.scatter(xlist, ylist, s =100)
ax2.streamplot(xp, yp, -np.gradient(p_mesh,axis=1),
-np.gradient(p_mesh,axis=0),
density=2)
for x,y,n in zip(xlist, ylist, names):
ax2.text(x+10, y+10, n)
ax2.scatter(xlist, ylist, s =100)
plt.show()
Уравнение производительности скважины можно записать в виде
\[ Q = T \Delta P J_D \tag{25}\]
где
\(Q\) - дебит жидкости скважины на поверхности, приведенный к стандартным условиям, м\(^3\)/сут. \[ Q = qB \]
\(T\) - параметр зависящий от гидропроводности пласта \[ T=\dfrac{18.41\mu B q }{\ k h} \tag{26}\]
\(\Delta P\) - депрессия на пласт, атм \[ \Delta P = \left(P_e-P_w\right) \tag{27}\]
\(J_D\) - безразмерный коэффициент продуктивности скважины, \[ J_D = \dfrac{1}{ \left(\ln{\dfrac{r_e}{r_w}} + S\right)} \tag{28}\]
Уравнение 25 можно интерпретировать следующим образом. Параметр \(T\) отвечает за свойства пласта и флюида на которые трудно повлиять в ходе эксплуатации. Это то, что дала природа в точке где находится скважина. Депрессия \(\Delta P\) – параметр которым можно управлять в ходе эксплуатации регулируя забойное давление. Например за счет установки насоса и задания параметров его работы. На этом параметре должно быть сосредоточено основное внимание при анализе работы скважины. Параметр \(J_D\) – определяет качество соединения скважины с пластом или качество заканчивания. Его мы можем выбирать при строительстве скважины и можем менять в ходе эксплуатации проводя ГТМ, хотя и достаточно большой ценой. Поскольку мы можем влиять на \(J_D\) важно понимать, какое оптимальное значение продуктивности можно достичь на конкретной скважине и как его можно изменить.
Задачей гидродинамических исследований является установление величин \(T\) и \(J_D\), хотя традиционно речь ведется об определении проницаемости \(k\) и скин-фактора \(S\).
Для совершенствования навыков работы с python выполните следующие задания:
Постройте график распределения давления в пласте для композитного пласта. В композитном пласте на расстоянии \(r<r_1\) проницаемость равна \(k=k_1\), а для \(r>=r_1\), \(k=k_2\).
Постройте двумерную тепловую или контурную карту распределения давления в пласте для моделей однородного и композитного пласта.
Рассчитайте среднюю величину давления в круговой области дренирования для однородного пласта. Насколько среднее давление в круговой области дренирования будет отличаться от давления на контуре. Чему будет равен коэффициент \(S\) в выражении \(Q=\dfrac{kh}{18.41\mu B} \dfrac{P_{res}-P_{wf}}{ln(\dfrac{r_e}{r_w})+S}\) при использовании вместо давления на контуре среднего давления? Постройте график, на котором будет отображаться распределение давления в зоне дренирования и величина среднего давления (в виде линии).
Для примера с несколькими скважинами имитирующими трещину ГРП рассчитайте дебиты скважин таким образом, чтобы забойное давление на всех скважинах было одинаковым. Постройте графики распределения давления в пласте. Постройте график дебитов вдоль “скважины”.
Для установившегося режима фильтрации давление в пласте не меняется. Для псевдо-установившегося режима постоянным остается перепад давления между пластом и забоем. После запуска, остановки или изменения режима работы скважины эти условия не выполняются. Давление в различных точках пласта может меняться по разному. Такой режим называют неустановившимся, а решения его описывающие нестационарными (зависят от времени).
Неустановившиеся решения уравнения фильтрации (transient solutions) представляют значительный практический интерес во многих задачах, включая задачи интерпретации ГДИС. В тоже время они относительно сложны и требуют применения компьютерных алгоритмов. В данном пособие проведение расчетов иллюстрируется с использованием python и макросов для Excel – Unifloc VBA.
Для решения уравнения фильтрации - линейного дифференциального уравнения в частных производных второго порядка необходимо задать начальные и граничные условия.
Самое простое решение можно получить для случая вертикальной скважины бесконечно малого радиуса запускающейся с постоянным дебитом. Условия соответствующие этому случаю можно выразить следующим образом:
В этом случае решение может быть выражено через функцию интегральной экспоненты \[ p_D(r_D,t_D) = - \frac{1}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right) \]
где \(-Ei(-x)\) - интегральная показательная функция, рисунок \(\ref{ris:Ei_plot_1}\).
\[ Ei(x)=-\int\limits_{x}^{\infty}\frac{e^{-t}}{t}\,\mathrm dt \]
# импортируем библиотеки для расчетов
# numpy используем для работы с массивами и подготовки данных для построения графиков.
# Также в некоторых функциях используем возможности векторных расчетов numpy
import numpy as np
# matplotlib используем для построения графиков
import matplotlib.pyplot as plt
# scipy.special используем как альтернативный вариант расчета специальных функций
import scipy.special as sc
# для скорости и удобства используем sc.expi
# построим ветки функции для положительных и отрицательных аргументов раздельно
x = np.arange(1e-5,3,0.01)
x1 = np.arange(-3,-1e-5,0.01)
plt.plot(x, sc.expi(x))
plt.plot(x1, sc.expi(x1))
plt.title("График функции Ei")
plt.xlabel("x")
plt.ylabel("Ei(x)")
plt.show()
Часто для проведения расчетов, особенно с использованием компьютерных библиотеке расчетов, бывает удобнее пользоваться модифицированной интегральной показательной функцией \(Ei_1(x)\) или \(E_1(x)\) или \(Ei_n(x)\) при \(n=1\). \[ Ei_n(x) = \int\limits_{1}^{\infty}\frac{e^{-tx}}{t^n}\,\mathrm dt \]
График интегральной показательной функции \(Ei_1(x)\) приведен на Рисунок 14. Для вещественных положительных \(x\in\mathbb R, x>0\) верно \(E_1(x) = - Ei( -x)\)
Функцию интегральной экспоненты можно представить в виде ряда.
\[ Ei(x)=-\int\limits_{x}^{\infty}\frac{e^{-t}}{t}\,\mathrm dt=\gamma+\operatorname{ln}|-x|+\sum\limits_{n\ge1}\frac{{-x}^n}{n!\cdot n}, \; x\in\mathbb R,\; \]
# зададим логарифмическое распределение точек вблизи нуля для построения графика
x = np.logspace(-10, 1, 100)
plt.rcParams["figure.figsize"] = (8,3)
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(x, sc.expi(-x))
ax1.plot(x, np.log(x) + 0.57721566481)
#ax1.set_title('Сравнение Ei и ln для обычных координат')
ax2.plot(x, sc.expi(-x))
ax2.plot(x, np.log(x) + 0.57721566481)
ax2.set_xscale('log')
#ax2.set_title('Сравнение Ei и ln для полулогарифмических координат')
plt.show()
Из приведенного выражения можно сделать выводы, что для маленьких значений аргумента функция интегральной экспоненты \(E_1(x)\) может быть аппроксимирована логарифмической зависимостью.
\[ E_1(x) = -ln(x) - \gamma \]
График сравнения функций \(E_1(x)\) и \(ln(x)\) показан на Рисунок 15. Видно, что хорошей аппроксимация будет только для маленьких значений аргумента \(x < 0.01\). Но для решения уравнения фильтрации именно эта зона представляет наибольший интерес.
Представление интегральной экспоненты в виде логарифмической аппроксимации удобно на практике, так как логарифм легче вычислять. В большинстве языков программирования и инструментов для проведения расчетов расчет логарифма реализован по умолчанию. А для расчета интегральной экспоненты, часто приходится предпринимать дополнительные шаги.
Решение уравнения фильтрации для линейного стока с учетом логарифмической аппроксимации можно представить в виде
\[ p_D(r_D,t_D) = \frac{1}{2} \left( ln \left( \dfrac{ t_D }{r_D^2} \right) +0.809 \right) \]
при использовании данного уравнения, следует помнить, что приближенное решение применимо при \(\dfrac{r_D^2}{4t_D} < 0.01\)
Решение линейного стока в размерных переменных
\[ p\left(r,t\right)=p_i-\frac{18.41q_sB\mu}{kh}\left(-\frac{1}{2}Ei\left(-\frac{\varphi\mu c_tr^2}{0.00144kt}\right)\right) \]
Решение с учетом логарифмической аппроксимации в размерных переменных
\[ p\left(r,t\right)=p_i-\frac{9.205q_sB\mu}{kh}\left(ln{\frac{kt}{\varphi\mu c_tr^2}}-7.12\right) \]
верно при \[ \frac{kt}{\varphi\mu c_tr^2}>70000 \]
Решения приведены для практических метрических единиц измерения, что можно увидеть по размерному коэффициенту.
Нестационарное решение с учетом скин-фактора будет иметь вид
\[ P(r, t) = P_{i} - \frac {9.205 {q_s} B\mu }{k h}(\ ln\frac {k t}{ \varphi \mu {c_t} {r^2}} -7.12 + 2S) \]
Уравнение фильтрации для радиального потока в линеаризованном виде можно записать в виде
\[ \frac{\partial p}{\partial t} = 0.00036 \dfrac{k}{\phi \mu c_t} \dfrac{1}{r} \left[ \dfrac{\partial}{\partial r} \left( r \dfrac{ \partial p} {\partial r} \right) \right] \]
Напомним, здесь
Часто для анализа уравнений неустановившейся фильтрации используются безразмерные переменные. Мы будем использовать переменные в виде:
\[ r_D = \frac{r}{r_w} \] \[ t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2} \] \[ p_D = \frac{kh}{ 18.41 q_s B \mu} \left( p_i - p \right) \]
Здесь использование единицы измерения СИ.
Использование безразмерных переменных позволяет упростить уравнение фильтрации, которое примет вид
\[ \frac{\partial p_D}{ \partial t_D} = \dfrac{1}{r_D} \left[ \dfrac{\partial}{\partial r_D} \left( r_D \dfrac{ \partial p_D} {\partial r_D} \right) \right] \]
Решение этого уравнения - функция безразмерного давления от безразмерных времени и расстояния $p_D(r_D, t_D) $
Для решения уравнения фильтрации - линейного дифференциального уравнения в частных производных второго порядка необходимо задать начальные и граничные условия. Самое простое решение можно получить для случая вертикальной скважины бесконечно малого радиуса запускающейся с постоянным дебитом. Условия соответствующие этому случаю можно выразить следующим образом
начальное условие. До запуска скважины в момент времени \(t_D = 0\) давление в пласте равно начальному во всех точках \(p=p_i\) \[ t_D < 0, p_D = 0 \]
условие постоянства дебита на скважине - граничное условие на скважине \[ \lim_{r_D \to 0} {r_D \frac{\partial p_D}{\partial r_D}} = -1 \]
условие на бесконечном расстоянии возмущения от скважине нет \[ r_D = \infty, p_D = 0 \]
В этом случае решение может быть выражено через функцию интегральной экспоненты \[ p_D(r_D,t_D) = - \frac{1}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right) \]
где -Ei(-x) - интегральная показательная функция.
Решение в размерных переменных можно записать как \[ p\left(r,t\right)=p_i-\frac{18.41q_sB\mu}{kh}\left(-\frac{1}{2} Ei \left(-\frac{\varphi\mu c_tr^2}{0.00144kt}\right)\right) \]
import scipy.special as sc
"""
Решение линейного стока уравнения фильтрации
"""
def pd_ei(td, rd):
"""
Решение линейного стока уравнения фильтрации
rd - безразмерное расстояние
td - безразмерное время
"""
return -1/2*sc.expi(-rd**2 / 4 / td)"""
построим графики решения линейного стока
в безразмерных переменных
"""
rd_arr = np.logspace(1, 3, 100)
td_arr = np.logspace(1, 3, 100)
# при построении используем векторный расчет
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,5))
for td in [1000, 10000, 100000]:
ax1.plot(rd_arr, pd_ei(td, rd_arr), label=f"td = {td}" )
ax1.set_title("Решение уравнения линейного стока для td = {}".format(td))
ax1.set_xlabel("rd")
ax1.set_ylabel("pd")
ax1.legend()
for rd in (1, 5 , 50):
ax2.plot(td_arr, pd_ei(td_arr, rd), label=f"rd = {rd}" )
ax2.set_title("Решение уравнения линейного стока для rd = {}".format(rd))
ax2.set_xlabel("td")
ax2.set_ylabel("pd")
ax2.legend()
plt.show()
Для расчетов в размерных переменных необходимо определить функции перевода размерных координаты в безразмерные и обратно. В целом расчеты станут более громоздкие.
"""
определим функции для перевода размерных переменных в безразмерные и обратно
пригодится для построения графиков и ведения расчетов
при наименовании функций придерживаемся следующих соглашений
сначала идет название того, что считаем
в конце указывается размерность результата, если это уместно
"""
def r_from_rd_m(rd, rw_m=0.1):
"""
перевод безразмерного расстояния в размерное
rd - безразмерное расстояние
rw_m - радиус скважины, м
"""
return rd*rw_m
def rd_from_r(r_m, rw_m=0.1):
"""
перевод размерного расстояния в безразмерное
r_m - размерное расстояние, м
rw_m - радиус скважины, м
"""
return r_m/rw_m
def t_from_td_hr(td, k_mD=10, phi=0.2, mu_cP=1, ct_1atm=1e-5, rw_m=0.1):
"""
перевод безразмерного времени в размерное, результат в часах
td - безразмерное время
k_mD - проницаемость пласта, мД
phi - пористость, доли единиц
mu_cP - динамическая вязкость флюида, сП
ct_1atm - общая сжимаемость, 1/атм
rw_m - радиус скважины, м
"""
return td * phi * mu_cP * ct_1atm * rw_m * rw_m / k_mD / 0.00036
def td_from_t(t_hr, k_mD=10, phi=0.2, mu_cP=1, ct_1atm=1e-5, rw_m=0.1):
"""
перевод размерного времени в безразмерное
t_hr - размерное время, час
k_mD - проницаемость пласта, мД
phi - пористость, доли единиц
mu_cP - динамическая вязкость флюида, сП
ct_1atm - общая сжимаемость, 1/атм
rw_m - радиус скважины, м
"""
return 0.00036 * t_hr * k_mD / (phi * mu_cP * ct_1atm * rw_m * rw_m)
def p_from_pd_atma(pd,
k_mD=10, h_m=10,
q_sm3day=20, b_m3m3=1.2, mu_cP=1, pi_atma=250):
"""
перевод безразмерного давления в размерное, результат
в абсолютных атмосферах
pd - безразмерное давление
k_mD - проницаемость пласта, мД
h_m - мощность пласта, м
q_sm3day - дебит на поверхности, м3/сут в с.у.
fvf_m3m3 - объемный коэффициент нефти, м3/м3
mu_cP - динамическая вязкость флюида, сП
pi_atma - начальное давление, абсолютные атм.
"""
return pi_atma - pd * 18.41 * q_sm3day * b_m3m3 * mu_cP / k_mD / h_m
def pd_from_p(p_atma,
k_mD=10, h_m=10,
q_sm3day=20, b_m3m3=1.2, mu_cP=1, pi_atma=250):
"""
перевод размерного давления в безразмерное
p_atma - давление
k_mD - проницаемость пласта, мД
h_m - мощность пласта, м
q_sm3day - дебит на поверхности, м3/сут в с.у.
fvf_m3m3 - объемный коэффициент нефти, м3/м3
mu_cP - динамическая вязкость флюида, сП
pi_atma - начальное давление, абсолютные атм.
"""
return (pi_atma - p_atma) / (18.41*q_sm3day*b_m3m3*mu_cP) * k_mD * h_m """
Построим графики распределения давления и изменения давления
в размерных координатах
"""
# исходные данные для построения графиков
h_m=10
q_sm3day=20
b_m3m3=1.2
mu_cP=1
pi_atma=250
rw_m=0.1
k_mD=10,
phi=0.2
ct_1atm=1e-5
r_arr = np.logspace(0.1, 3, 100)
rd_arr = rd_from_r(r_arr, rw_m=rw_m)
t_arr = np.linspace(1, 100, 100)
td_arr = td_from_t(t_arr,
k_mD=k_mD,
phi=phi,
mu_cP=mu_cP,
ct_1atm=ct_1atm,
rw_m=rw_m)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,5))
for td in td_from_t(np.array([1, 10, 50]),
k_mD=k_mD,
phi=phi,
mu_cP=mu_cP,
ct_1atm=ct_1atm,
rw_m=rw_m):
ax1.plot(r_arr,
p_from_pd_atma(pd=pd_ei(td,rd_arr),
h_m=h_m,
q_sm3day=q_sm3day,
b_m3m3=b_m3m3,
mu_cP=mu_cP,
pi_atma=pi_atma),
label=f"t = {t_from_td_hr(td,
k_mD=k_mD,
phi=phi,
mu_cP=mu_cP,
ct_1atm=ct_1atm,
rw_m=rw_m)[0]:.2f} час")
ax1.set_title("Воронка депрессии")
ax1.set_xlabel("r, m")
ax1.set_ylabel("p, atm")
ax1.legend()
for rd in rd_from_r(np.array([rw_m, 10, 50])):
ax2.plot(t_arr,
p_from_pd_atma(pd=pd_ei(td_arr, rd),
h_m=h_m,
q_sm3day=q_sm3day,
b_m3m3=b_m3m3,
mu_cP=mu_cP,
pi_atma=pi_atma),
label=f"r = {r_from_rd_m(rd, rw_m=rw_m)}" )
ax2.set_title("Изменение давления на расстоянии")
ax2.set_xlabel("t, час")
ax2.set_ylabel("p, atm")
ax2.legend()
plt.show()
Иногда для манипуляции со сложными выражениями удобно бывает использовать sympy
import sympy as sp
# объявим символьные переменные
r_d, p_d, t_d, gamma = sp.symbols('r_d p_d t_d gamma')
r, p, t, q, b, mu, phi, ct, k, h, p_i, s, r_w = sp.symbols('r p t q b mu phi ct k h p_i s r_w')
# зададим исходное уравнение
eq1 = sp.Eq(p_d, -1/2*sp.Ei(-r_d**2/4 / t_d))
eq1\(\displaystyle p_{d} = - 0.5 \operatorname{Ei}{\left(- \frac{r_{d}^{2}}{4 t_{d}} \right)}\)
решение с интегральной экспонентой можно задать используя scipy.special или сгенерировать из решения sympy
pd_ei_ = sp.lambdify([t_d, r_d],eq1.rhs, modules=["scipy","numpy"])
# сравним решения заданные явно и сгенерированные из решения `sympy`
print(pd_ei(1000,1))
print(pd_ei_(1000,1))
print(pd_ei(1000,1)==pd_ei_(1000,1))3.8585419797881815
3.8585419797881815
Truefig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,5))
rd_arr = np.logspace(0,2,100)
for td in [1, 10, 100, 1000, 10000]:
ax1.plot(rd_arr, pd_ei_(td, rd_arr))
ax1.plot(rd_arr, pd_ei(td, rd_arr))
rd_arr = np.logspace(0,2,100)
for td in [1, 10, 100, 1000, 10000]:
ax2.plot(rd_arr, pd_ei_(td, rd_arr))
ax2.plot(rd_arr, pd_ei(td, rd_arr))
ax2.set_xscale('log')
Связь размерных и безразмерных переменных задается в виде
\[ r_D = \frac{r}{r_w} \] \[ t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2} \] \[ p_D = \frac{kh}{ 18.41 q_s B \mu} \left( p_i - p \right) \]
Здесь использование единицы измерения СИ.
# подставим вместо безразмерных переменных размерные
eq2 = eq1.subs(r_d, r/r_w)
eq2 = eq2.subs(t_d, 0.00036*k*t/(phi*mu*ct*r_w**2))
eq2 = eq2.subs(p_d, k*h/(18.41*q*b*mu)*(p_i-p))
# решим полученное уравнение относительно перепада давления
eq3 = sp.solve(eq2,p_i-p)
# выведем первое решение
eq3[0]\(\displaystyle - \frac{9.205 b \mu q \operatorname{Ei}{\left(- \frac{694.444444444444 ct \mu \phi r^{2}}{k t} \right)}}{h k}\)
Решение в размерных переменных можно записать как \[ p\left(r,t\right)=p_i-\frac{18.41q_sB\mu}{kh}\left(-\frac{1}{2} Ei \left(-\frac{\varphi\mu c_tr^2}{0.00144kt}\right)\right) \]
# сгенерируем решение
pd_ei_full_ = sp.lambdify([t, r, k, h, mu, ct, phi, q, b, p_i],p_i-eq3[0],
modules=["scipy","numpy"])
# нарисуем воронку депрессии используя сгенерированное решение
r_arr = np.logspace(-1,2,100)
for tt in [0.1, 1, 10]:
plt.plot(r_arr,
pd_ei_full_(tt, r_arr,
k=10, h=10, mu=1, ct=1e-4, phi=0.2, q=10, b=1, p_i=250))
Решение с интегральной экспонентой может быть заменено приблеженным решением с использованием логарифма
\[ p_D(r_D,t_D) = - \frac{1}{2} \ln \left( \dfrac{ r_D^2}{4t_d} \right) - \frac{1}{2}\gamma \]
где \(\gamma = 0.57721566481\) - константа Эйлера
на графике от времени в полулогарифмических координатах логарифмическое приближение выглядит как кривая с наклоном \(0.5\)
упростим последнее выражение с использованием sympy
# объявим символьные переменные
r_d, p_d, t_d, gamma = sp.symbols('r_d p_d t_d gamma')
r, p, t, q, b, mu, phi, ct, k, h, p_i, s, r_w = sp.symbols('r p t q b mu phi ct k h p_i s r_w')
# зададим исходное уравнение
eq = sp.Eq(p_d, -1/2*sp.ln(r_d**2/4 / t_d) -1/2*gamma)
print('Исходное уравнение - логарифмическое приближение решения линейного стока')
display(eq)Исходное уравнение - логарифмическое приближение решения линейного стока\(\displaystyle p_{d} = - 0.5 \gamma - 0.5 \log{\left(\frac{r_{d}^{2}}{4 t_{d}} \right)}\)
# подставим значение для gamma
eq = sp.simplify(sp.expand(eq.subs(gamma, 0.57721566)))
eq\(\displaystyle p_{d} = 0.404539350559945 - 0.5 \log{\left(\frac{r_{d}^{2}}{t_{d}} \right)}\)
pd_ei_ln_ = sp.lambdify([t_d, r_d],eq.rhs, modules=["scipy","numpy"]) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,5))
rd_arr = np.logspace(0,2,100)
for td_ in [1, 10, 100, 1000, 10000]:
ax1.plot(rd_arr, pd_ei_(td_, rd_arr))
ax1.plot(rd_arr, pd_ei_ln_(td_, rd_arr))
rd_arr = np.logspace(0,2,100)
for td_ in [1, 10, 100, 1000, 10000]:
ax2.plot(rd_arr, pd_ei_(td_, rd_arr))
ax2.plot(rd_arr, pd_ei_ln_(td_, rd_arr))
ax2.set_xscale('log')
# подставим вместо безразмерных переменных размерные
eq = eq.subs(r_d, r/r_w)
eq = eq.subs(t_d, 0.00036*k*t/(phi*mu*ct*r_w**2))
eq = eq.subs(p_d, k*h/(18.41*q*b*mu)*(p_i-p))
# решим полученное уравнение относительно перепада давления
eq1 = sp.solve(eq,p_i-p)
# выведем первое решение
eq1[0]\(\displaystyle \frac{9.205 \cdot 10^{-14} b \mu q \left(80907870111989.0 - 100000000000000.0 \log{\left(\frac{2777.77777777778 ct \mu \phi r^{2}}{k t} \right)}\right)}{h k}\)
упростим решение
eq5 = sp.simplify(sp.expand(eq1[0]))
display(eq5)\(\displaystyle \frac{b \mu q \left(- 9.205 \log{\left(\frac{ct \mu \phi r^{2}}{k t} \right)} - 65.5426176327539\right)}{h k}\)
Для совершенствования навыков работы с python выполните следующие задания:
Один из самых простых примеров применения суперпозиции. Предполагаем, что добывающая скважина в однородном изотропном пласте запускается в момент времени t=0 и работает t_p_hr часов, после чего останавливается. После остановки скважины забойное давление растет - и мы получим кривую восстановления давления.
Пусть решение задачи запуска скважины (падения давления) будет \(P_D(t_D, r_D)\). Тогда решение для изменения давления при запуске и последующей остановки скважины можно представить в виде \[\begin{equation} P_{bu.D}(t_D, t_{prod.D}, r_D) = P_D(t_D) - P_D(t_D-t_{prod.D}, r_D) \cdot \mathcal{H}(t_D-t_{prod.D}) \end{equation}\]
где
\(t_D\) - безразмерное время после запуска скважины,
\(t_{prod.D}\) - безразмерное время работы скважины после запуска
\(\mathcal{H}\) - ступенчатая функция Хевисайда (в некоторых книгах обозначается как \(\theta\))
\(P_D(t_D, r_D)\) - безразмерное давление - решение задачи запуска скважины (падения давления)
\(P_{bu.D}(t_D, t_{prod.D}, r_D)\) - безразмерное давление- решение задачи запуска скважины и последующей остановки скважины
Для проведения векторных расчетов в python удобно выражение с использованием функции Хевисайда
\[ \mathcal{H} = \begin{cases}0 & x < 0\\1 & x = 0\\1 & x > 0\end{cases}\]
Применение функции Хевисайда позволяет избежать в расчетных функциях применение условных операторов в явном виде для отдельных элементов входных массивов. Это потенциально ускоряет расчет.
def pd_build_up(td, td_p, rd):
"""
расчет давления для запуска и последующей остановки скважины
td - время после запуска
td_p - время безразмерное - которое скважина работала до остановки
rd - расстояния от скважины
"""
# применение функции Хевисайда здесь делает расчет корректным
# для входных векторов td
return pd_ei(td, rd) - np.heaviside(td-td_p,1) * pd_ei(td-td_p, rd)t_arr = np.logspace(-10, 2, 1000)
t_prod_hr = 24
k = 10 # проницаемость
q = 30 # дебит
# переведем размерный массив времени в безразмерные величины
# некоторые параметры можно было бы пропустить и оставить значения по умолчанию
# но для полноты приведем их в явном виде
td_arr = td_from_t(t_arr,
k_mD=k, phi=0.2, mu_cP=1, ct_1atm=1e-05, rw_m=0.1)
td_prod = td_from_t(t_prod_hr,
k_mD=k, phi=0.2, mu_cP=1, ct_1atm=1e-05, rw_m=0.1)
print('время работы скважины {:.2f} часа, что соответсвует безразмерному времени {:.2f}'.format(t_prod_hr, td_prod))
# для заданного массива безразмерных времен рассчитаем безразмерные давления
pd_arr = pd_build_up(td_arr, td_prod, rd=1)
# построение графика
plt.rcParams["figure.figsize"] = (8,3)
plt.plot(td_arr, pd_arr)
plt.xlabel('td')
plt.ylabel('pd')
plt.title('Решение в безразмерных координатах')
plt.show()время работы скважины 24.00 часа, что соответсвует безразмерному времени 4320000.00
# переведем безразмерные координаты в размерные
p_arr = p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1, pi_atma=250)
plt.rcParams["figure.figsize"] = (8,3)
plt.plot(t_arr, p_arr)
plt.xlabel('t, hr')
plt.ylabel('p, atma')
plt.title('Решение в размерных координатах')
plt.show()
r_arr = np.logspace(1, 3, 100)
t_arr = np.logspace(-1, 2, 10)
t_arr[0] = 0.001
t_prod = 100
tv, rv = np.meshgrid(td_from_t(t_arr),rd_from_r(r_arr))
pd_arr =pd_ei(tv, rd=rv)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = [12,5])
fig.suptitle('Изменение давления в пласте при запуске скважины')
ax1.plot(r_arr,
p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax1.set_label(t_arr)
ax1.set_xlabel('r, m')
ax1.set_ylabel('$p_{wf}, atma$')
ax2.plot(r_arr,
p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax2.set_xscale('log')
ax2.set_xlabel('r, m')
ax2.set_ylabel('$p_{wf}, atma$')
plt.show()
pd_arr_bu =pd_ei(tv + td_from_t(t_prod), rd=rv) - pd_ei(tv, rd=rv)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = [12,5])
fig.suptitle('Изменение давления в пласте при остановке скважины')
ax1.plot(r_arr,
p_from_pd_atma(pd_arr_bu,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax1.set_label(t_arr)
ax1.set_xlabel('r, m')
ax1.set_ylabel('$p_{wf}, atma$')
ax2.plot(r_arr,
p_from_pd_atma(pd_arr_bu,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax2.set_xscale('log')
ax2.set_xlabel('r, m')
ax2.set_ylabel('$p_{wf}, atma$')
plt.show()
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = [12,5])
fig.suptitle('Изменение давления в пласте при запуске и остановке скважины')
ax1.plot(r_arr,
p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax1.plot(r_arr,
p_from_pd_atma(pd_arr_bu,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax1.set_label(t_arr)
ax1.set_xlabel('r, m')
ax1.set_ylabel('$p_{wf}, atma$')
ax2.plot(r_arr,
p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax2.plot(r_arr,
p_from_pd_atma(pd_arr_bu,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax2.set_xscale('log')
ax2.set_xlabel('r, m')
ax2.set_ylabel('$p_{wf}, atma$')
plt.show()
Для расчета изменения давления при переменном дебите введем произвольное референсное значение дебита $ q_{ref} $ (например первое не нулевое значение дебита при запуске скважины). Используем это значение для определения безразмерного давления. \[ p_D = \frac{kh}{ 18.41 q_{ref} B \mu} \left( p_i - p \right) \]
и безразмерного дебита
\[ q_D = \frac{q}{q_{ref}} \]
Тогда, используя принцип суперпозиции, можем выписать выражение для изменения давления на скважине и вокруг нее для произвольного момента времени
\[ P_{mr.D}(t_D, r_D) = \sum_i \left[ q_{D(i)}-q_{D(i-1)} \right] \cdot p_D\left(t_D-t_{D(i)}, r_D\right)\cdot \mathcal{H}(t_D-t_{D(i)}) \]
где
# создадим историю изменения дебитов
t_history = np.array([ 0., 2., 24. ], dtype=np.float64)
q_history = np.array([10., 5., 0.], dtype=np.float64)
# массивы должны быть одной и той же длины
def pd_superposition(td, td_hist, qd_hist, rd):
"""
расчет безразмерного давления для последовательности
безразмерных дебитов
td - время расчета после запуска, безразмерное
td_hist - массив времен изменения режимов работы скважин, безразмерное
qd_hist - массив дебитов установленных после изменения режима работы,
безразмерное
"""
# принудительно добавим нули во входные массивы,
# чтобы учесть запуск скважины
qdh = np.hstack([0, qd_hist])
tdh = np.hstack([0, td_hist])
# построим дебиты виртуальных скважин -
# разности реальных дебитов при переключении
delta_qd = np.hstack([0, np.diff(qdh)])
# векторная магия - время как вектор и переключения дебитов тоже
# организуем сумму по временам, внутри сумма по переключениям
# делаем при помощи расчета meshgrid и поиска накопленных сумм
qd_v, td_v =np.meshgrid(delta_qd, td)
# используем куммулятивную сумму numpy для того что суммировать
dpd = np.cumsum(qd_v * pd_ei((td_v - tdh),
rd) * np.heaviside((td_v - tdh), 1), 1)
# последний столбец - полная сумма, которая нужна в качестве результата
return dpd[:,-1]
def q_superposition(t, t_hist, q_hist):
"""
расчет давления для запуска и последующей остановки скважины
t_hr - время после запуска в часах
t_hist_hr - массив времен изменения режимов работы скважин
q_hist_sm3day - массив дебитов установленных после изменения режима работы
"""
def interpolate_constant(x, xp, yp):
indices = np.searchsorted(xp, x, side='right')
y = np.concatenate(([0], yp))
return y[indices]
q=[]
for ti in t:
q.append(interpolate_constant(ti, t_hist, q_hist))
return qtd_arr = np.linspace(1e-3, 9e4, 1000)
td_history = np.array([ 0., 2e4, 5e4 ], dtype=np.float64)
qd_history = np.array([10., 5., 0.], dtype=np.float64)
plt.rcParams["figure.figsize"] = (12,5)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(td_arr,
q_superposition(td_arr, td_history, qd_history))
ax2.plot(td_arr,
pd_superposition(td_arr, td_history, qd_history, rd=1),
color='red')
ax2.set_xlabel('td')
ax1.set_ylabel('qd')
ax2.set_ylabel('pd')
plt.show()
def p_superposition_atma(t_hr, t_hist_hr, q_hist_sm3day,
k_mD=10, h_m=10, b_m3m3=1.2, mu_cP=1, pi_atma=250,
phi=0.2, ct_1atm=1e-05, rw_m=0.1):
"""
расчет давления для запуска и последующей остановки скважины
t_hr - время после запуска в часах
t_hist_hr - массив времен изменения режимов работы скважин
q_hist_sm3day - массив дебитов установленных после изменения режима работы
k_mD=10 - проницаемость, мД,
h_m=10 - мощность пласта, м,
b_m3m3=1.2 - объемный коэффициент, м3/м3,
mu_cP=1 - вязкость нефти, сП,
pi_atma=250 - начальное давление, атм,
phi=0.2 - пористость, доли единиц,
ct_1atm=1e-05 - общая сжимаемость, 1/атм,
rw_m=0.1 - радиус скважины
"""
q_ref=1.
td = td_from_t(t_hr, k_mD=k_mD, phi=phi, mu_cP=mu_cP,
ct_1atm=ct_1atm, rw_m=rw_m)
td_hist = td_from_t(t_hist_hr, k_mD=k_mD, phi=phi, mu_cP=mu_cP,
ct_1atm=ct_1atm, rw_m=rw_m)
return p_from_pd_atma(pd_superposition(td, td_hist, q_hist_sm3day / q_ref, 1),
k_mD=10, h_m=10,
q_sm3day=q_ref, b_m3m3=1.2, mu_cP=1, pi_atma=250)t_arr = np.arange(1e-6, 50, 1e-2)
plt.rcParams["figure.figsize"] = (12,5)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(t_arr, q_superposition(t_arr, t_history, q_history))
ax2.plot(t_arr, p_superposition_atma(t_arr, t_history, q_history), color='red')
ax2.set_xlabel('t, hr')
ax1.set_ylabel('q, m3/day')
ax2.set_ylabel('p, atma')
plt.show()
t_history = np.array([0, 2, 10, 24, 24.1,24.2,24.3,24.4,24.5 ])
q_history = np.array([10, 11, 12, 11, 10, 9, 8, 7, 6])
plt.rcParams["figure.figsize"] = (12,5)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(t_arr, q_superposition(t_arr, t_history, q_history))
ax2.plot(t_arr, p_superposition_atma(t_arr, t_history, q_history, 1), color='red')
ax2.set_xlabel('t, hr')
ax1.set_ylabel('q, m3/day')
ax2.set_ylabel('p, atma')
plt.show()
rd, td, qd = sp.symbols('r_D t_D q_D')
# запишем решение в символьном виде для постоянного дебита
eq = -qd/2 * sp.Ei(-rd **2 / 4 / td)
display(eq)\(\displaystyle - \frac{q_{D} \operatorname{Ei}{\left(- \frac{r_{D}^{2}}{4 t_{D}} \right)}}{2}\)
для линейно меняющегося дебита \(dq_D t_D\) решение можно представить как интеграл
\[ p_D = \int_0^{t_D}{- \frac{dq_D}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right) dt_D} \]
# проинтегрируем решение по времени
eq1 = sp.integrate(eq, td)
display(eq1)\(\displaystyle - \frac{q_{D} \left(\frac{r_{D}^{2} \operatorname{Ei}{\left(- \frac{r_{D}^{2}}{4 t_{D}} \right)}}{4} + t_{D} \operatorname{Ei}{\left(- \frac{r_{D}^{2}}{4 t_{D}} \right)} + t_{D} e^{- \frac{r_{D}^{2}}{4 t_{D}}}\right)}{2}\)
Для линейно меняющегося дебита во времени (как и для любой другой зависимости) надо решение проинтегрировать по времени (надо бы подробнее расписать - сделать это позже, например как у Щелкачева в основах нестационарной фильтрации на стр 321).
Для линейной зависимости дебита от времени \[ Q_D = dQ_D \cdot t_D \] можно получить выражение \[ p_D(r_D,t_D, dQ_D) =-\frac{dQ_D t_D }{2} \left[ \left( 1+ \frac{r_D^2}{4 t_D} \right) Ei \left(- \dfrac{r_D^2}{4t_D} \right) + e^{-\dfrac{r_D^2}{4t_D}} \right] \]
где \(dQ_D\) - скорость изменения дебита.
Для таблично заданных дебитов и времен можно оценить
\[ dQ_{D(i)} = \dfrac{Q_{D(i)}-Q_{D(i-1)}}{t_{D(i)} - t_{D(i-1)} } \]
Cравните формулу (11) с формулой (9.68) в книге Щелкачева “Основы неустановившейся фильтрации”
Тогда, используя принцип суперпозиции, можем выписать выражение для изменения давления на скважине и вокруг нее для произвольного момента времени
\[ P_{mr.D}(t_D, r_D) = \sum_i p_D\left(t_D-t_{D(i)}, r_D, dQ_{D(i+1)} - dQ_{D(i)}\right)\cdot \mathcal{H}(t_D-t_{D(i)}) \]
где
следует обратить внимание, при суперпозиции скорость изменения дебита вычисляется как \(dQ_{D(i+1)} - dQ_{D(i)}\). при реализации расчета необходимо предусмотреть, чтобы для первого и последнего шага расчет прошел корректно. Для этого можно, например, добавить к массивам дебитов и времени дополнительный значения в начале и в конце массивов соответствующие постоянным значениям дебита.
Также надо учитывать, что в приведенном выражении массивы должны начинаться со значений \(Q_D=0\)
# Решение линейного стока уравнения фильтрации
def pd_lin_ei(td, rd=1, dqd_dtd=1):
"""
Решение линейного стока уравнения фильтрации
rd - безразмерное расстояние
td - безразмерное время
"""
# при расчете убедимся, что td=0 не повлияет на расчет,
# даже если td массив и нулевой только один элемент
td = np.array(td, dtype = float)
pd = (1 + rd**2/4/td) * (-sc.expi(-rd**2/4/td)) - np.exp(-rd**2/4/td)
return dqd_dtd * td * pd / 2
def pd_superposition_lin(td, td_hist, qd_hist):
"""
расчет безразмерного давления для последовательности
безразмерных дебитов
- td - время расчета после запуска, безразмерное
- td_hist - массив времен изменения режимов работы, безразмерное
- qd_hist - массив дебитов после изменения режима работы, безразмерное
"""
# принудительно добавим нули во входные массивы,
# чтобы учесть запуск скважины
qdh = np.hstack([qd_hist])
tdh = np.hstack([td_hist])
# построим дебиты виртуальных скважин
# разности реальных дебитов при переключении
delta_qd = np.hstack([np.diff(qdh),0])
delta_td = np.hstack([np.diff(tdh),1])
dq_dt = delta_qd / delta_td
dq_dt = np.diff(np.hstack([0, delta_qd / delta_td]))
# референсный безразмерный дебит это 1
# векторная магия - время как вектор и переключения дебитов тоже
# организуем сумму по временам, внутри сумма по переключениям
# делаем при помощи расчета meshgrid и поиска накопленных сумм
qd_v, td_v =np.meshgrid(delta_qd, td)
dpd = np.cumsum(pd_lin_ei((td_v - tdh),
dqd_dtd=dq_dt) * np.heaviside((td_v - tdh), 1),1 )
return dpd[:,-1]
def p_superposition_lin_atma(t_hr, t_hist_hr, q_hist_sm3day,
k_mD=10, h_m=10, b_m3m3=1.2, mu_cP=1, pi_atma=250,
phi=0.2, ct_1atm=1e-05, rw_m=0.1):
"""
расчет давления для запуска и последующей остановки скважины
t_hr - время после запуска в часах
t_hist_hr - массив времен изменения режимов работы скважин
q_hist_sm3day - массив дебитов установленных после изменения режима работы
k_mD=10 - проницаемость, мД,
h_m=10 - мощность пласта, м,
b_m3m3=1.2 - объемный коэффициент, м3/м3,
mu_cP=1 - вязкость нефти, сП,
pi_atma=250 - начальное давление, атм,
phi=0.2 - пористость, доли единиц,
ct_1atm=1e-05 - общая сжимаемость, 1/атм,
rw_m=0.1 - радиус скважины
"""
q_ref = 1.
td = td_from_t(t_hr,
k_mD=k_mD, phi=phi, mu_cP=mu_cP, ct_1atm=ct_1atm, rw_m=rw_m)
td_hist = td_from_t(t_hist_hr,
k_mD=k_mD, phi=phi, mu_cP=mu_cP,
ct_1atm=ct_1atm, rw_m=rw_m)
return p_from_pd_atma(pd_superposition_lin(td, td_hist, q_hist_sm3day / q_ref),
k_mD=10, h_m=10,
q_sm3day=q_ref, b_m3m3=1.2, mu_cP=1, pi_atma=250)td_arr = np.linspace(1e-3, 2000, 2000)
td_history = np.array([0., 1, 50, 100, 300, 500 ], dtype=np.float64)
qd_history = np.array([0., 3, 1.4, 1.0, 1.01, 1], dtype=np.float64)
plt.rcParams["figure.figsize"] = (12,8)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(td_arr, q_superposition(td_arr, td_history, qd_history))
ax1.plot(td_arr, np.interp(td_arr, td_history, qd_history))
ax2.plot(td_arr,
pd_superposition_lin(td_arr, td_history, qd_history),
color='red')
ax2.plot(td_arr,
pd_superposition(td_arr, td_history, qd_history, rd=1),
color='blue')
ax2.set_xlabel('td')
ax1.set_ylabel('qd')
ax2.set_ylabel('pd')
plt.show()
t_history = np.array([0, .02, 10, 24, 24.1,24.2,24.3,24.4,24.5 ])
q_history = np.array([0, 11, 12, 12, 10, 9, 8, 7, 6])
plt.rcParams["figure.figsize"] = (12,5)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(t_arr, q_superposition(t_arr, t_history, q_history))
ax1.plot(t_arr, np.interp(t_arr, t_history, q_history))
ax2.plot(t_arr,
p_superposition_atma(t_arr, t_history, q_history),
color='red')
ax2.plot(t_arr,
p_superposition_lin_atma(t_arr, t_history, q_history),
color='green')
ax2.set_xlabel('t, hr')
ax1.set_ylabel('q, m3/day')
ax2.set_ylabel('p, atma')
plt.show()
Нестационарное решение в безразмерных переменных \[ p_D(r_D,t_D) = - \frac{1}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right) \] где безразмерные переменные введены как \[ r_D = \frac{r}{r_w} \] \[ t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2} \] \[ p_D = \frac{kh}{ 18.41 q_s B \mu} \left( p_i - p \right) \]
Здесь использование единицы измерения СИ. - \(r_w\) - радиус скважины, м - \(r\) - расстояние от центра скважины до точки в пласте, м - \(q_s\) - дебит скважины на поверхности, приведенный к нормальным условиям м3/сут - \(\phi\) - пористость, доли единиц - \(\mu\) - вязкость нефти в пласте, сП - \(B\) - объемный коэффициент нефти, м3/м3 - \(p_i\) - начальное давление в пласте, атм - \(p\) - давление на расстоянии \(r\), атм - \(c_t\) - общая сжимаемость системы в пласте, 1/атм
Для этих же безразмерных переменных, считая начальное давление равным давлению на контуре можно записать стационарное решение для движения в круговом пласте
\[ p_D = \ln r_{eD} - \ln r_D \]
сравним это решение с логарифмической аппроксимацией (1)
\[ p_D(r_D,t_D) = - \frac{q_D}{2} \left[ \ln \left( \dfrac{ r_D^2}{4t_d} \right) +\gamma \right] \]
которое можно преобразовать к виду \[ p_D(r_D,t_D) = - q_D \ln r_D + \frac{q_D}{2} \left[ \ln(4t_D) -\gamma \right] \]
сравнивая со стационарным решением можно найти выражение безразмерного радиуса контура в зависимости от безразмерного времени \[ \ln r_{eD} = \frac{1}{2}(\ln(4t_D)-\gamma) \]
\[ r_{eD} = \sqrt { 4t_D e^{-\gamma} } \]
# оценим значение величины под корнем
print(4*np.exp(-0.57721566481))2.2458379344731085наконец получим \[ r_{eD} = \sqrt {2.2458 t_D} \]
это значение называют радиусом влияния скважины. Используя это значение для определенного момента времени можно получить стационарное распределение давления в системе хорошо приближающее решение линейного стока работающего в бесконечном пласте. Можно считать это расстояние за расстояние на которое распространяется влияние скважины.
достижение радиуса влияния внешних границ будет обуславливать начало перехода от неустановившегося режима фильтрации к режиму обусловленному влиянием границ - стационарному для границы постоянного давления или псевдоустановившемуся для границы непротекания.
# Решение линейного стока уравнения фильтрации
def pd_ei(td, rd):
"""
Решение линейного стока уравнения фильтрации
- rd - безразмерное расстояние
- td - безразмерное время
"""
return -1/2*sc.expi(-rd**2 / 4 / td)
def pd_ss(rd, red):
"""
стационарное решение в безразмерных переменных
- rd - безразмерное расстояние
- red - безразмерное расстояние до границы
"""
return np.log(red/rd)# зададим точки расстояний для отрисовки графика
rdl = np.logspace(-1, 3 , 100)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = [12,5])
# построим первый график в обычных координатах
for td in [10, 100, 1000]:
ax1.plot(rdl, pd_ei(td, rdl))
red = (td*2.2458)**0.5
ax1.plot(rdl, pd_ss(rdl, red))
# построим второй график в полулогарифмических координатах
for td in [10, 100, 1000]:
ax2.plot(rdl, pd_ei(td, rdl))
red = (td*2.2458)**0.5
ax2.plot(rdl, pd_ss(rdl, red))
plt.xscale('log')
plt.show()
Преобразования размерных величин удобно выполнять с модулем символьных вычислений python - sympy. Преобразования размерностей ключевых величин полезно знать наизусть, хотя всегда можно найти их в таблицах. Значения многих физические константы зашины в модуле scipy.constants, иногда это оказывается удобным, при этом автоматически будет учитываться достаточно большое количество знаков после запятой в константах. Рассмотрим размерности ряда величин широко применяемых в нефтяном инжиниринге.
# импортируем sympy и константы из scipy
import sympy as sp
import scipy.constants as const
# в модуле scipy.constants есть значения общепринятых констант -например значение pi
print(const.pi)3.141592653589793В СИ измеряется в [м\(^3\)/сек], в практических метрических единицах измеряется в [м\(^3\)/сут], в американских промысловых единицах измеряется в [bbl/day].
# выведем некоторые переводные коэффициенты для объемных расходов
print(f'Одни [сут] = {24*60*60} = {const.day} [сек]')
print(f'Один [м3/сут] = {1/const.day} [м3/сек]')
print(f'Один баррель в день [bbl/day] = {const.bbl} [м3/сут]')
print(f'Один баррель в день [bbl/day] = {const.bbl/const.day} [м3/сек]')
print(f'Один [м3/сут] = {1/const.day} [м3/сек]')Одни [сут] = 86400 = 86400.0 [сек]
Один [м3/сут] = 1.1574074074074073e-05 [м3/сек]
Один баррель в день [bbl/day] = 0.15898729492799998 [м3/сут]
Один баррель в день [bbl/day] = 1.8401307283333331e-06 [м3/сек]
Один [м3/сут] = 1.1574074074074073e-05 [м3/сек]В СИ измеряется в [м\(^2\)], в практических метрических единицах измеряется в [мД], в американских промысловых единицах измеряется в [mD].
Определение: в пористой среде с проницаемостью в один Дарси для поддержания течения жидкости с динамической вязкостью 1 сП со скоростью фильтрации 1 см/с необходимо поддерживать перепад давления жидкости приблизительно в одну атмосферу на 1 см вдоль направления течения. При использовании физической атмосферы для расчета перепада давления (физическая атмосфера = 101 325 Па) единица проницаемости равняется приблизительно 0.986923 мкм².
В отечественной литературе при определении дарси в качестве величины атмосферы было принято использовать техническую атмосферу (1 кгс/см² = 98 066,5 Па), так что для величины дарси получалось значение приблизительно 1,02 мкм², причём эпизодические случаи использования западного определения дарси специально отмечались [ru.wikipedia.org/wiki/Дарси]. Согласно ГОСТ 26450.2-85 величины 1 Дарси \(= 0.9869⋅10^{−12}\) м².
print(f'Один [мД] = {1e5/const.atm * 1e-15} [м²]')Один [мД] = 9.86923266716013e-16 [м²]AT = 98066.5 # technical atmosphere in Pa, техническая атмосфера в Па
print(f'Один [psi] в [Па] = {const.psi}')
print(f'Один [bar] в [Па] = {const.bar}')
print(f'Один [atm] в [Па] = {const.atm}')
print(f'Один [at] в [Па] = {AT}')
print(f'Один [atm] в [psi] = {const.atm/const.psi}')Один [psi] в [Па] = 6894.757293168361
Один [bar] в [Па] = 100000.0
Один [atm] в [Па] = 101325.0
Один [at] в [Па] = 98066.5
Один [atm] в [psi] = 14.69594877551345Используя рассчитанные выше переводные коэффициенты для различных размерных величин рассчитаем переводной коэффициент в формуле Дюпюи
\[ Q = \frac{ 2 \pi kh}{ \mu B} \frac{ \left( p_i - p \right) } {\ln{\dfrac{r_e}{r_w}} +S } \]
# зададим переменные sympy
Q, k, h, mu, B, pres, pwf, re, rw, S, pi = sp.symbols('Q k h mu B p_res p_wf r_e r_w S pi',
real=True,
positive=True)
# определим уравнение
eq = sp.Eq(Q, 2 * pi * k * h / (mu * B) * (pres - pwf) / (sp.ln(re/rw) + S))
display(eq)\(\displaystyle Q = \frac{2 h k \pi \left(p_{res} - p_{wf}\right)}{B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
eq = eq.subs(Q, 1/const.day * Q) # дебит, [м3/сут] в [м3/сек]
f_k = 1e5/const.atm * 1e-15
eq = eq.subs(k, f_k * k) # проницаемость, [мД] в [м2]
eq = eq.subs(mu, 1e-3 * mu) # вязкость, [сП] в [Па сек]
eq = eq.subs(pres, const.atm * pres) # давление [атм] в [Па]
eq = eq.subs(pwf, const.atm * pwf) # давление [атм] в [Па]
eq = eq.subs(pi, const.pi)
display(eq)\(\displaystyle 1.15740740740741 \cdot 10^{-5} Q = \frac{6.20102176874373 \cdot 10^{-12} h k \left(101325.0 p_{res} - 101325.0 p_{wf}\right)}{B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
Решим полученное уравнение относительно Q и упростим средствами sympy
eq1 = sp.simplify(sp.solve(eq,Q)[0])
display(sp.Eq(Q,eq1))\(\displaystyle Q = \frac{0.0542867210540316 h k \left(p_{res} - p_{wf}\right)}{B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
Выделим полученную константу в явном виде и найдем обратную величину - это и будет необходимый нам переводной коэффициент.
f = 1/eq1.args[0]
print(f)18.4207110060064По умолчанию sympy автоматически организует порядок элементов в своих выражениях. Этот порядок может отличаться от привычного - хотя и суть формул при этом не меняется. Применяя некоторые хитрости можно заставить sympy вывести выражения в приемлимом виде.
a = sp.symbols('a')
eq2 = eq1.subs(eq1.args[0],1/a)
with sp.evaluate(False):
display(eq2.subs(a, f))\(\displaystyle \frac{h k \left(p_{res} - p_{wf}\right)}{18.4207110060064 B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
Но иногда результат проще переписать руками в нужном виде. В итоге уравнение Дюпюи в практических метрических единицах измерения примет вид.
\[ Q = \frac{ kh}{18.42 \mu B} \frac{ \left( p_i - p \right) } {\ln{\dfrac{r_e}{r_w}} +S } \]
где